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ABSTRACT 


A procedure is described to determine lateral buckling 
loads for initially straight beams. Loading and beam geom- 
etry may vary along the length. Warping rigidity is not 
considered. End conditions of considerable generality may 
be treated. The algorithm depends on solving a convergent 
sequence of sixteenth order elgenproblems. A computer pro- 
gram implementing this procedure has been developed and is 


described herein. 
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Meee bbe nonee LELCATIONS 


The coordinate system used is a right-handed orthogonal 
X,y,z system with positive senses as indicated in figure l. 
The longitudinal centroidal axis of the beam lies along the 
+z direction. The principal axes of a normal cross section 
are Originally in the x and y directions but after torsional 
displacement (rotation) they are in the € and n directions, 


Sieetigure 2 (page 10). 





Figure 1. General configuration and coordinate system. 


The beam displacements are defined by three variables: 
vertical deflection v, OTS WORSE. deflection u, and angular 
rotation ». Positive senses for v and u are as indicated in 
figure 1. The positive sense for ¢ is indicated in figure 2. 
The deformations u, v and $ are considered to be "small" in 


a mathematical sense, so that, for example sin ¢ w ¢. 








/ 


/ 


Pigure 2. Positive sense of $¢ and source of bending 
moment N-oM. 

Defining the problem requires specifications for five 
physical characteristics of the beam. They are beam length 
Meet lexural rigidity for bending in the vertical plane STL. 
Meerued! rigidity for bending in the horizontal plane EI2, 
memeiendal rigidity C, and warping rigidity Cl. The last 
four may vary along the length of the beam. 

Torsional rigidity C is defined by the following equa- 


en 


Ce OleOkROUE) / (ANGLE OF TWIST PER UNIT LENGTH) (1) 


Powe simplicity if Metation, we use numerical suffixes 
rather than numerical subscripts. 


M6 





Warping rigidity Cl is related to nonuniform torsion. A 
discussion here of this property would be too digressive; 
See section 5.3 of reference (2). 

It 1s the underlying purpose of this entire analysis to 
determine the total externally applied load P which will 
Brice incipient lateral buckling. The only loading con- 
Sidered is distributed load in the negative y direction. 
igestype of loading Cuniform, triangular or whatever) is 


mpeeitied beforehand by a dimensionless function f(z) such 


that 
wi Se Gay eee (2) 


where the load shape function w(z) has dimensions of force 
per unit length. Concentrated loads and concentrated moments 
can be approximated by appropriate large local variations in 
imegoe)6fhe function f(z) will be denoted by x5 when it 
appears later. The constant R has dimensions of force times 
distance squared. The magnitude of R is arbitrary and is 
chosen at the convenience of the user. 

The load function has been defined. We are interested 
mmeetermining the multiplying factor Q such that lateral 
ferme does mot take place if the actual loading is less 
Piemmenw 2) and does take place if the actual loading exceeds 
Qw(z). The entire purpose of the program developed from 
this analysis is to determine this multiplier Q. In order 


to obtain nondimensionalized equations, it 1S convenient to 


de 





Mimeoduce FP the total applied load for incipient buckling, 


L 
P = a | wz) dz (3) 
0 


defined by 


per that on PL{/R (5) 


For the special case where the load contributions in the 
negative y direction are canceled by the load contributions 
ie qe positive y direction the following alternate form of 
equation 3 is applied. 

L 


P “| Oncz) | dz (6) 


0 


The load may be applied above, on, or below the cen- 
troidal axis as specified by the load eccentricity, e(2), 
See cigure 6 (page 20). 

For notational and programming convenience a series of 
dimensionless functions x are defined. The first six are 
used to specify beam properties and loading. A seventh is 


Simply a function having constant unit value. 


ale Z 





cl. = ELL/R Ce) 


= ER (8) 
Soa nc R (9) 
x4 = RL“/Cl (10) 
x5 = w(z) L°/R = f(z) (11) 
x6 = e(Z)/L 2) 
ao eeaal les) 


The nondimensionalizing scheme used to define the x 
Functions is applied throughout the problem development. 
The dimensionless forms for vertical and horizontal deflec- 


tion are 
rl ee, © 1 Gl) 
and | 32S hye eS.) 


where 8 1S a dimensionless function which ultimately will be 
constructed from several functions. 


We will use the dimensionless parameter 
G Seite Gia) 


to specify axial position. Accordingly we will think of the 
ferumnetions as functions of ¢, thus x5 = x5(c). Differentia- 
tion with respect to z and f are indicated as (= = Y') and 


(Ss = M2) respectively. The two operations are related by 


the equation 


Gig?) 


Wes 





The developments which follow require considerable inte- 
gration. Two notational devices are used to facilitate both 
the readability and the transcription of equations. Integra- 


ifemom 1S Slenified by p, thus 


px3 | ecules yal(e 


Pyemuation of any function at the right end of the beam, i.e., 
for 5 = 1, is indicated by the prefix *. Several examples 


illustrating this notation are 


1 
ape | ale ale (18) 
0 
Bx5x9 = x5(1) * x9(1) (19) 
XL 
x2 px8 = 2a | x8(t)dz (20) 
0 
1 G 
tox5x6px1l = f asconocert [ sarcorasras (21) 
0 0 


(In the last equation ® is a dummy variable.) 
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w(z) w(z) 


rc Bq [pL 


k3v (0) V(z) k4v(L) 
Figure 3. Diagram of forces and moments in the vertical 
plane. 


Our analysis is restricted to beams which are uncon- 
strained except at one or both ends. There are six possible 
constraints at each end each of which is modeled as a spring. 
The spring constant k describing such a constraint can 
actually vary from zero (perfectly free, no constraint at 
all) to infinity (perfect or ideal constraint) inclusive. 
Therefore to avoid mathematical difficulties: the analysis 
employs dimensionless spring parameters, designated by the 
symbol a (with an appropriate identifying suffix). Fora 
torsion spring, such as that corresponding to the moment 
klv'(0) shown in figure 3 (page 15), the relations between 
me) and al are 


al Kiel aR /L) CoS, 


kek ely er (ilea = ced) (23) 


Note that al = 0 corresponds to kl = 0 (perfect freedom) 


and al = 1 corresponds to kl = ~ (perfect constraint). 


iS 





Similarly, for a linear spring such as spring number 3, the 


mMelations are 


ad = Ratee 2 Rhu) (24) 


ene a Cll = aa) (25) 


Figure 5 (page 13) shows that the point of application of a 
linear spring is not limited to the centroidal axis. The 
vertical deflection constraints may be applied at distances 
cl and c2 above the centroidal axis, and the horizontal de- 
flection constraints may be applied at distances c3 and c4 
above the centroidal axis. The positive sense is as indi- 
cated in figure 5, and the dimensionless forms of these dis- 


tances are 


Gl = el/L Cs) 
G20= e727 L 27) 
Gan = e37 ks (28) 
G4 = c4/L (29) 


The physical and geometrical parameters used in the prob- 
lem have now been defined. They will be used in the struc- 


tural analysis developed in the next section. 


les 


Peo nochuURAL ANALYSIS 


There are three structural deformations examined in the 
analysis. They are vertical bending, horizontal bending, 
and torsion. 

The vertical bending problem is disjoint from the remain- 
der of the analysis, and it 1s solved first. The notation 
used in defining vertical bending will not be shown because 
the problem is simple and familiar. It is sufficient to say 
that a set of four simultaneous nonhomogeneous equations are 
formed, representing the constraint conditions at each end 
which pertain to bending in the vertical (yz) plane. This 
system is then solved by a standard algebraic analysis, and 
the result 1s used to create three dimensionless functions 


for shear, moment, and deflection. 


Shear Ome = V7 © (30) 
Moment x9 = M/PL Ca) 
Beptectton xO = v/OL G2) 


The remainder of the analysis involves studying the com- 
Bined horizontal bending and torsion problem. Figure 4 
(page 18) shows the positive sense for the constant horizon- 
tal shear H, the linearly varying horizontal bending moment 


ie eand horizontal displacement u. 


Le 





k5u' (0 - c 5 of [te 


k7Lu(0)+c19(0)] J|>° 


Figure 4. Diagram showing forees and moments 
Mieenewnorrazontal plane, 


The appropriate dimensionless variables are 


Bees wll) eo: Ob Ee Se) 
Bh = H/P ; H = PB4 = BYP (34) 
n = N/PL; N= PLn (35) 


Dimensionless shear B4 is one of ultimately eight un- 
Pm@fewm scalar quantities (Bl, B2,...., B8) which will be in- 
troduced in the process of determining Q. The alternate 
forms indicated in the preceding equations are intended to 
familiarize the reader with notational forms which will be 
freely employed in later developments. 

The illustration of the rotated cross section in figure 


2 (page 10) shows that the moment causing bending about the 


nN axis is N cos @ - M sin $6, which reduces to N - Md by small 
angle approximation. From elementary bending theory we have 
EI2 u" = N - @M (36 ) 
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and by substitution from equation 7 it may be written as 


B = Qx2(n - ox9) 


37) 


Examination of the simple statics of figure 4 (page 18) 


yields 
Zee COD. = NCZ) 
In dimensionless form this becomes 
pier temeo = 1 Co) 


and may be rewritten as 












n = BYpx7 + BSx7 
n€Q) = B5 
k4v(L) 
k3v(0) 
k5u' (0) k6u'(L) 


M(Q) 
ktv' (0) 


k7Lu(0)+c16(0) J 


| 
v(Q) 
k9 (0) 


ise 
Cro 


M(L) 
k2v'(L) 


so) te 


CBR, 


3) 


(40 ) 


C41) 


KaCu(L)+c20(L)] 


c TCL) 


mesure 5, Force and moment diagrams of end sections. 


i 





Figures 5 and 6 (page 20) show the forces and moments 


used to develop the torsion moment equation 


T(0)-H[v(z)-v(0)] + VOz)uCz)-V(0)u(0) 


T(z) = 
Z 
‘| MEeZDLUGzZ te (zo (Zz) |dz Gu?) 
0 
Introducing the dimensionless function 
(43) 


tls. 1/ PL 


and using additional dimensionless functions previously de- 


equation 31 may be written as 





mee d , 
f= oO) —OBY [xl0—x10(0) ]+x86-x8(0)8(0)+pxS(Bt+x6¢o) (44) 
V(z),down 
V(0),u0 2? 
7 
b. CENTRAL PORTION 
———ew | (0) WL) ee 
= me 99 (0) K106 (CL) mnie 
GZ amet (06360) VILIC4O (La ess 
a. LEFT END c. RIGHT END 
mesure 6, Diagram of left, central, and right 


portions for torsion analysis. 


ZC 





Reference (1) develops the following third order equa- 


[meme which accounts for the warping rigidity of the beam. 
Cie = eo. + T - u'M (45) 
In dimensionless form this may be written as 
@ = x4¥(Qt - Qx9B + x36) (46) 
The function x4 has previously been defined as 
x4 = RL“/Cl (47) 


Cl is the coefficient of the highest (i.e., third) order 
derivative in equation 45, and if Cl becomes small, very 
troublesome mathematical difficulties are encountered, cf. 
Reference (3). Accordingly, the analysis for the case Cl 
= 0 cannot be obtained from that for Cl # 0, given in refer- 
ence (1), except with the greatest of difficulty. It is 
for this reason that the separate analysis represented by 
this thesis, was undertaken. The analysis herein, to this 
point, is; as has been said before, essentially identical 
to that in reference (1). From this point on there are 
essential differences. 

The present analysis presumes that Cl = 0, so that equa- 


iemom 35 takes the form 
‘ aa : 
Gr Oxce werdp = t) (48) 


The inverse of x3, Oe oceurs frequently in what 
follows. The term is clumsy in this form and will be re- 


designated as 
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Y = x3 (49) 
Therefore equation 48 is written as 
6 = QOY(x9B - t) (50) 


Two of the end constraint conditions of the original 
problem, namely those relating to warping constraints at the 
end, are not available for and do not appear in the present 
specialized problem. 

The original problem had six constraints at each end. 
The present problem has five constraints at each end. Four 
constraint conditions were employed in solving the problem 
Semeending in the vertical plane. Thus a total of six con- 
straint conditions remain at this stage of the development. 
The first four result from the constraints against horizon- 
tal deflection and end rotations in the x-z plane. The last 
two result from constraints against rotation about the z 


axis. 


a5B(0) = QCl-a5) n(0) = 0 (51) 
ae) + OC l-a6) n(1) = 0 (52) 
wae +Glo(d))] + OBY(l=-a7) = 0 -(53) 
Gite 2o(1) i) = OBe(i-as) = 0 (54) 
a96(0)+Q(1-a9)[t(0)+63x8(0)6(0)-B4YG1] = 0 (55) 
Po y= OC l=ald){t(1)+G4x8(1)¢(1)-BuG2] = 0 (56) 


Although this notation appears formidable, it is con- 


venient, explicit, and unequivocal. The terms are 
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sufficiently complicated that no matter what notation might 
Pema carerul attention to detail is required. 

There 1s no need to demonstrate the complete development 
of the six boundary equations, but for illustrative purposes 
the derivations of two typical constraint equations, 51 and 
Soeewill be given. 

The constraint against rotation of the left end about a 
vertical (y) axis, as shown in figure 6 (page 20) results 


in the following equation. 
oueGO ye NGO) =. 0 a7) 


Employing the dimensionless variable 8 = u/L and on = N/PL, 
Meme equation 17, writing the spring constant k5 in terms 


of the spring parameter a5, viz, 
oe = mao Gl—as ) (58) 
and recalling equation 5, it 1s easy to obtain the equation 
-a58(0)-Q(1-a5)n(0) = 0 (59) 


ims 1S equation 51. 

Equation 53 defines the constraint against horizontal 
deflection of the left end of the beam. The deflection is 
u(0) + c1o(0), see figure 5 (page 19), where the second term 
results from permitting the spring to be attached at a dis- 
tance cl above the centroidal axis. The constraint equation 
is 


vom ciao 7] + ft = 0 (60) 


ZnS 





where as previously defined H = PB4Y, cl = LGl, and 
k7 = a7R/L°(1-a7) om) 


Using the dimensionless quantities previously introduced, 
it 1S a Simple matter to put equation 60 in the form ex- 
hibited as equation 53. 

The structural analysis has now been completed and is re- 
Meesenced by equations 37, 44, 50, 51, 52, 53, 54, 55, and 
56. The mathematical method of solution will be addressed 


mex. 
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IV. MATHEMATICAL ANALYSIS 


Before going into the details of the analysis, we will 
try to describe it in general terms. 

An expression is formed for b which involves an unknown 
mumeeron (xii in what follows) the initial form of which is 


arbitrarily assumed except for the normalizations 


i 
© 


oil NG Ole. Caz 


‘pl xl1| Ces) 


il 
| 


The corresponding function $$ 1s obtained by integration and 
equation 37 is used to get 8. This 1s followed by two more 
integrations to get 8 and 8. These are employed in the tor- 
Seem equation 44 to get the function t. Finally this is 

used in equation 48 to get a new expression for . Each in- 
tegration introduces a new unknown scalar quantity. In one 
Way or another a total of eight such unknown scalar constants 
Bl, B2,...,B8 are introduced. The six constraint equations, 
plus two more which will be introduced in what follows, pro- 
vide a set of eight linear homogeneous equations in these un- 
known B's which thus leads to an eigenproblem of eighth 
order. The parameter Q appears in the coefficients to the 
meest and second powers. The quadratic eighth order eigen- 
problem is transformed to a linear sixteenth order eigenprob- 
lem. Since, generally, both zero and infinite elgenvalues 


are contained in the solution, a special, relatively new 
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algorithm, the QZ method, capable of dealing with such sys- 
tems, 1S employed for the solution of this eigenproblem. 
The appropriate eigenvalue is chosen and the corresponding 
eigenvector provides the solution for the unknown B's so 
that the new xll can be uniquely determined. This new xll 
is used in place of the original xll and the process 1s re- 
peated until there is convergence. No proof is offered 
that indeed convergence must ultimately take place; because 
of the complexity of the problem such a proof might be very 
difficult to establish. However, experience with the pro- 
cedure described here indicates that convergence may indeed 
be expected. The details of this procedure are now given in 
what follows. 


The function ¢ is initially specified as 
@ = Blxll + B2x7 (64) 


where xll is an arbitrarily assumed function satisfying equa- 
femems 62 and 63 and where Bl and B2 are the third and fourth 
of eight unknown constants that will be employed. The first 
was B4 from equation 34 and the second was BS from equation 
41. 


Integrating equation 64 yields 
(= remap? to Box (65) 


and when equations 40 and 65 are substituted into equation 37 


we arrive at 
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8 = -B1Qx2x9pxll - B2Qx2x9px7 - B30x2x9 


+BUQx2px7 + BSQx2 (66) 
Integrating twice more, we get 


8 = -BlQpx2x9pxll - B2Qpx2x9px7 - B3Qpx2x9 


weyODs 2px? + BaOpX2 + BEQX7 (67) 
fen > LOppx2x3pxll = B2Qppx2x3px7 - B3Qppx2x9 
moO pex2p xe soCppx2 + B6Qpx?7 + B7Qx7 (68) 


For convenience the next step is to rewrite equation 44 


as 
t = a8 + px5x6o - QBYX - Qx8(0)B7 + QBS (69) 
where Te OU) 2 Oley (70) 
SS sel) Seta) Gia) 
dee ape o to 8 (72>) 


The symbol q denotes an operator acting on the function 8, 
Meveea function which is multiplied by 8. 

When the preceding equations, evaluated at 7@=0 or f=1 as 
the case may be, are substituted in boundary Bone cene) oma 
tions 51 through 56, a system of six linear simultaneous 
homogeneous equations in eight unknown constants B is obtained. 

In the previous section, which addressed the structural 
analysis, equations 51 and 53 were derived for illustrative 
purposes. Now for illustrative continuity these two equa- 
tions will be used to show the method of formulation of the 


Six linear simultaneous homogeneous equations. 


a 





Recalling that equation 51 is 
a5B(0) - Q(l-a5) n(0) = 0 
and substituting 


8(0) = QB6 CNS 


moe 


Bo (74) 
we now factor out Q leaving 
~B5(l-a5) + B6aS = 0 (75) 
The second illustrative equation, 53, is 
-a7[B(0) + G1¢(0)] + QOBYCl-a7) = 0 


and upon substituting 


CD Ke (76) 
Bqu) = B70 OA) 

we arrive at 
Boclay Spaewoel—ajy) + B/Qa7 = 0 Ci3)) 


If equations 67 and 69 are substituted into equation 50 


Sieerollowing "new" » is created. 
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> = ¥[B1(-0°x9px2x9px1l i 0° appx2x9pxl1l =. Oe pede joy clea, 


as 


+ 


. 


B2 (-0° x9px2x9px7 + Q¢ appx2x9px7 - Qpx5x6px7) 
B3(-Q°x9px2x9 23 Q* appx2x9 - Qpxs5x6 ) 

Bu (Q°x9px2px7 i Q°x - 0° appx2px7) 

BS (0° x9px2 - Q- appx?) 

B6 (Q°x9 - Q° apx7) 


Bel OF YI (79) 


(The coefficient of B7 vanishes identically.) 


To complete the eilgensystem involving the eight unknown 


B's, two additional equations are needed. 


The seventh. equation, consistent with equation 62, is 


new ¢(0) - B2 = 0 (80) 


The eighth equation is 


(*DO) O14 =~ C0 ae = 0 Goa 


and is arrived at by observing that if 6 given in equation 64 


and created from an assumed xll, were, somehow, a "correct" 


' then it and the $ recovered in equation 79 would be iden- 


tically equal to each other. The integration in equation 81 


achieves a certain "smoothing" or averaging which has been 


found to be preferable to simply equating the new and old 


functions at an arbitrary value of ¢. 


The coefficients of the elgensystem equations generally 


myo lLve 


Q to the zeroth, first, and second powers so that 


the system may be exhibited in the form 


(D + EQ + FQ*)B = 0 (82) 


ZS 





The nonzero elements of matrices D, E, and F are 


meee aos DiG = a5, D21 = -a6*px2x9pxl1l, 
Waee= =do*pxz2xdpx’, D23 = -a6*px2x9, D24 = a6“px2px7-a6tl1, 


D25 = a6“px2-a6t+t1, D26 = a6, D332 Gla7, 


D4l = a8G2*pxll, D42 = a8G2, D43 = a8sG2, 


Dece= a9, D6l = al0*pxll, D62 = al0, D63 = ald, 


Pee 1, DSi = -"*pxll, D82 = -1, E34 = 1-a7, 
Be? = a/, EY] = -a8*ppx2x9px 11, 
E42 = -a8"ppx2x9px7, E43 = -a8*ppx2x9, 


E44 = a8*ppx2px7+a8-1, EYS = a8*ppx2, 

Becr= a8, EN7 = a8, E53 = (1-a9)G3x8(0), 

E54 = (a9-1)Gl, E6l = (al0-1)[*px5x6pxl1+G4*x8(1)pxll], 
Boe] (al0=-1)[*px5x6px7 + G4*x8(1)], | 

E63 = (al0-1)[*px5x6+G4*x8(1)], E64 = (1-al0)G2, 

meee" DYpxoxOpxll, E82 = -*“pYpxox6px/7, 

Bes = —*pYpxox6, F58 = l-ag, 

moe (l—al0)*qppx2x9pxll, F62 = (1l-al0)*qppx2x9px7, 
lemma 1 —-al10)*qppx2x9, FEY = (1-al0)(*X-*qppx2px7), 


momma l0-l)*qppx2, F66 = (Cal0-1)*qpx7, 


meee = al0=1, F776 = x9(0)Y(0), F78 = -Y(0), 
meme —"pYxopx2xdpxll + “pYqppx2x9pxl1l, 
P82 = =*pYx9px2x9px7 + *pYqppx2x3px7, 

Moe e= —="pDYX9DX2x9 + “DYqppx2x?3, 

F84 = *pYx9px2px7 + *pYX - *pYqppx2px7, 
feos) *>DYx9px2 — “pYqppx2 , 

moor “pYxd = “pYqpx/ , 

F88 = -*pY 
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[For reference purposes all the preceding equations bear the 
monmber 84.) 
This eighth order quadratic eigensystem may be expanded 


to a sixteenth order linear eigensystem in the form 


VZ = QWZ (85) 
E 
where Ves (86) 
to 8 
-F QO 
= Coe) 
Oe del 
QB 
Ze 4 (88) 


and I and 0 respectively represent the eighth order unit and 
null matrices. 

The elgenproblem has now been defined in a linear form, 
but because the matrices V and W may be singular not just 
any method will satisfactorily solve the problem. We have 
chosen the QZ method described in reference (3). From the 
Sixteen eigenvalues produced the smallest positive eigenvalue 
is selected. The last eight elements of the associated eigen- 
vector are the desired constants Bl through B8. 

Using equation 79 a new 6 is now formed and from it a 
new xll is created by means of the following normalization 


process. 


x11 = [4 - CO I/L*p | > - ¢ 60) | (89) 


new new QUS new he 


Sak 


The iteration process is now repeated until satisfactory 
convergence is attained. The converged smallest positive 
eigenvalue Q is now used to determine the buckling load from 


either equation 3 or equation 6. 
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ecw len eer DEMENTATION OF THEORY 


To this point the functions x have been treated as con- 
tinuous functions for use in an analytic solution. However, 
it is obvious that solving the problem, even for simple 
m@memwmons, by analytic means 1S prohibitively difficult. 

To overcome this problem a numerical method of solution was 
employed. The beam was divided into a finite number of 
equal sections. The functions x were then redefined as vec- 
mem@~emwith Fortran designations x(1,I), x(2,I)..., ete. 

fmemme £F = 1, 2, 3,...-,N. The elements of any such vector 
represent function values at the N points of subdivision of 
the beam into equal subsections of length L/(N-1). The 
vector x(10,I) was used twice, first to represent x10 and 
then to store X of equation 71. 

Beeevectors in this form readily lend themselves to 
manipulation. The simple subroutines developed to perform 
the operations are, with the exception of the integration 
Sememe, mot worthy of note here. 

Trapezoidal integration was first used during the debug- 
ging of the raw, untested program. This method gives satis- 
Mmaeeory results with an error of anes where h is the in- 
terval length. Once it was established that the program 
would give satisfactory solutions the refinement process was 
begun. A higher order numerical integration scheme was em- 


ployed next. It is described by Milne in reference (5), 
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and results in an error Gey 

The QZ subroutine is a canned program in the library of 
the Naval Postgraduate School Computer Center. It was trans- 
scribed from reference (6), with very little change, by 
Mr. Roger Hilleary. To date, it has been the only program 
that has given satisfactory solutions to the present eigen- 
value problem. 

The complete theory described in the preceding sections 
of this thesis has been implemented in the form of a sub- 
routine LATBRO, a complete listing of which is contained in 
Appendix A. 

The user must write a main program which supplies dimen- 
Sloning statements, number of intervals, specifications for 
vectors xl through x6, specifications for the constants G, 
and the ten spring parameters. The following is an example 
of such a main program. It was used to solve what is listed 
as Case 6, Table 1, page 37. 

IMPLICIT REAL*8(A-H,0-Z) 
heer Cea Ol) -ALFACI2),G(4),S,ONE, ZERO 
MIREGER KPC10) 
COMMON X,ALFA,G,N,KP 
N=101 
ZLERO=0.D+0 
ONE=1.D+0 
Dod f=1.N 
X(1,I)=ONE 
G2. 1 )=ONE 
X(3,1I)=ONE 
X(5,1)=ZERO 
X(6,I)=ZERO 
econ TINUE 
%(5,1)=ONE 
DO 2 I=1,4 


pe) = ZERO 
DO 3 I=1,12 
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3 ALFA(I)=ZERO 
ALFA(2 )=ONE 
ALFAC4)=ONE 
ALFA(6 )=ONE 
ALFA(8)=ONE 
ALFA(10)=ONE 


KE 3S 

CALL LATBRO(S) 
SLOP 

END 


crs 





ieee Ceol aNGeyAND VERIFICATION 


As mentioned in the previous section, the numerical in- 
tegration methods that have been employed have their short- 
comings. 

If the integrand is "smooth}' either the trapezoidal or 
the higher order Milne method may be used. The latter 1s 
preferable because of its greater accuracy. 

Both methods introduce errors, however, when there are 
steps, impulses, or doublets to be integrated. We have 
found that the Milne and trapezoidal methods can be used to 
solve problems with a concentrated load. The error depends 
linearly on the reciprocal of the number of elements into 
which the beam is divided and extrapolation may be used to 
obtain an excellent result. As yet we do not have experi- 
ence with concentrated moment loadings (doublets). This 
subject is discussed in more detail in reference (1). 

The following table shows the cases that have been tes- 
ted to date. In each of these cases all the G's and the 
eccentricity x6 were zero. Each of these has been success- 
ell 

The values for Q shown in the last column of the table 
Peeeemombtained with x2 = 1, x3 = 1, and by simple scaling 


analysis lead to total load 


P= 0 Varney ake (83 ) 
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for lateral buckling. The values reported agree with values 
to be found in reference (2). 

The following are details of several cases. 

Case l 

AS shown, it is a cantilever beam fully restrained at the 


right end. The loading is uniform. 





Figure 7. Illustration of side view of Case 1 beam. 


The subroutine produced a converged value of Q of 12.854. 
Timoshenko's result is 12.85; cf. page 261 ref (2). 

Case 2 

This is simply Case 1 turned end for end. It was used 
to test additional equations and also resulted in a converged 
Seated! load output of 12.854 which is identical to the re- 


eee for Case 1, as it should be. 
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ease 5 
The converged critical load output was 16.936 and Timo- 


shenko's result is 16.94; ef. page 269 ref (2). 


Figure 8. Top and side views Case 5 beam. 


In several of the cases listed in the table above, the 
ieee does not interact with one or more of the springs, 
so that the corresponding spring parameters could have an 


arbitrary value without changing the results. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


The subroutine LATBRO in its present form should be able 
to solve beam problems with beams that have uniform or vari- 
able cross sections, with loading on the centroidal axis, 
and with spring constraints attached at the centroidal axis. 
Loading that involves concentrated moments, as mentioned be- 
fore, has not yet been tested. 

The program was designed to solve problems where loading 
and spring constraints may act away from the centroidal axis. 


It 1s recommended that follow on testing begin with 


i loading off the centroidal axis. 

2. spring constraints attached away from the centroidal 
axis. 

Ss» concentrated moment loading. 


It is expected, with the exception of minor programming 
errors, that problems with any combination of the above will 


be readily solved by the subroutine. 
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Listing of LATBRO and Adjunct Subroutines 
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SOLVE THE X,Z BENDING PROBLEM 
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THE FIRST STEP WILL BE TO NORMALIZE THE LOADING VECTGR X5. 
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CALL NORV(5) 
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NOW PREPARE THE COEFFICIENT ARRAYS FOR THE X,Z BENDING PROBLEM. 
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NOW USING A CANNED SUBROUTINE TO SOLVE THe LINEAR SYSTEM. 
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NOW WORKING ON TERMS DIVIDED BY X3 BUT NOT INVOLVING X11. 
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